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Abstract. We present a computational study for a family of discontinuous Galerkin meth- 
ods for the one dimensional Vlasov-Poisson system, recently introduced in g]. We introduce 
a slight modification of the methods to allow for feasible computations while preserving the 
properties of the original methods. We study numerically the verification of the theoretical 
and convergence analysis, discussing also the conservation properties of the schemes. The 
methods are validated through their application to some of the benchmarks in the simulation 
of plasma physics. 



Numerical simulation has become a major tool for understanding the complex behavior of 
a plasma or a particle beam in many situations. This is due not only to the large number of 
physical applications and technological implications of the behavior of plasmas, but also to 
the intrinsic difficulties of the models used to describe such behavior. In fact, it was recog- 
nized long time ago that there does not exist any fully satisfactory macroscopic model (fluid 
equations) which can be used to describe the particle interaction in laser-fusion problems. 
In contrast, microscopic models (kinetic equations) can provide a more accurate description 
of the plasmas. 

One of the simplest model problems that is currently used in the simulation of plasmas 
is the Vlasov-Poisson system. Such system describes the evolution of a plasma of charged 
particles (electrons and ions) under the effects of the transport and self-consistent electric 
fleld. The unknown, typically denoted by f{x^v^t) (with x standing for position, v for 
velocity and t for time), represents the distribution function of particles (ions, electrons, 
etc.) in the phase space. The coupling with a self-consistent electrostatic fleld (neglecting 
magnetic effects) is taken into account through the Poisson equation. The nonlinear structure 
of the system prevents from obtaining analytical solutions, except for a few academic cases 
(see the surveys ^35l [l3l [26] for a good description on the state of the art of the mathematical 
analysis of the problem). Therefore, numerical simulations have to be performed to study 
realistic physical phenomena. 

At the present time, there can be distinguished two main classes of numerical methods for 
simulating plasmas; Lagrangian (or probabilistic) and Eulerian (or deterministic) methods. 
The former class include all different types of particle methods [121 ESJ |47l [34[ [HI [32l U\ and 
has been a preferred choice since the beginnings of numerical simulations in plasma physics 
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in the 60's, due to their simphcity and low computational cost. The basic idea behind these 
methods is to approximate the motion of the plasma by a finite number of macro-particles 
in the phase space whose trajectories are computed from the characteristics of the Vlasov 
equation, while the electrostatic field is computed by collecting the charge density on a fixed 
mesh of the physical space. Although this class of methods represents a feasible option and 
potentially might allow for resolving the whole 3 + 3 + 1 dimensional problem, their inherent 
numerical noise precludes from obtaining an accurate description of the distribution function 
in the phase space in many interesting cases. This lack of precision can be overcome by using 
a method from the second class; an Eulerian solver. These type of methods are nothing but 
classical (or new) numerical schemes discretizing the Vlasov equation on a (fixed) mesh of 
the phase space. Among them, the most widely used are finite volumes [271 EH] and semi- 
lagrangian methods [3ll [IHl [I5l [HJ [23l EH [25l [9]. Finite volumes (FV) are a simple and 
inexpensive option, but in general are low order if one wants to retain the basic conservative 
properties of the scheme. 

Semi-lagrangian schemes (sometimes consider in-between Eulerian and Lagrangian solvers) 
have become a popular option, since they can achieve high order allowing at the same time 
for time integration with large time steps. However, special care in needed to compute the 
origin of the characteristics with high order interpolation without spoiling the local character 
of the reconstruction. A nice numerical study comparing some of the different methods use 
in plasma simulations is presented in [28] . 

In this paper we present a computational study with discontinuous Galerkin (DG) meth- 
ods. DG methods are finite element methods that use discontinuous polynomials. Their local 
construction endow the methods with good local conservation properties without sacrificing 
the order of accuracy. This is one of the main motivations for their use in plasma simula- 
tions. But it also provides the methods a built-in parallelism which allows for parallelization 
of the resulting algorithms. The methods have also many other attractive features: they are 
extremely fiexible in handling /ip-adaptivity, the boundary conditions are imposed weakly 
and the DG mass matrices are block-diagonal which results in very efficient time-stepping 
algorithms in the context of time-dependent problems, as it is the case here. In spite that 
nowdays, DG methods are consider for approximating problems of almost any kind, their 
use for kinetic equations, and more particularly for simulation of plasmas has only been 
contemplated very recently. Among the computational works, we mention the use of DG in 
a multi-waterbag approximation of the VP system in [TT]; a piecewise constant DG solver 
for VP in |33j (which require extremely fine meshes) and semi-lagrangian schemes combined 
with high order DG interpolation are presented in [45l|43]. In both works, the authors also 
use the positivity preserving limiter introduced in [48]. 

A theoretical work has been presented in [H [5] , where the authors have introduced and 
analyzed a family of semi-discrete DG schemes for the VP system with periodic boundary 
conditions, for the one and multi-dimensional cases, respectively. The authors show optimal 
error estimates for both the distribution function and the electrostatic field, and they study 
the conservation properties of the proposed schemes. Due to the local construction of the 
DG schemes, total mass (or charge) conservation is shown to hold easily. This property 
is essential in the numerical approximation to VP, since it is required for guaranteeing the 
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well-posedness of the related Poisson problem. The authors also propose a novel DG scheme 
that is shown to preserve the total energy of the VP system. Their proof however requires 
the assumption that the DG finite element spaces contain at least all quadratic polynomials. 

In this work, we undertake the issues of verification and validation of these family of DG 
schemes, for the one-dimensional VP system. To accomplish both tasks, we first discuss 
how the schemes can be efficiently implemented in practice, even in parallel. For the space 
discretization, the methods introduced in [4j are based on the coupling of a DG discretization 
for the Vlasov equation (transport equation) together with a mixed finite element (possibly 
discontinuous) approximation to the Poisson problem. Here, however, we present two slight 
variations of the DG approximation for the Vlasov equation, to allow for feasible compu- 
tations. The modifications are done in the definition of the numerical fiux involving the 
coupling with the approximate electrostatic field (hence the nonlinearity) . The definition 
in [4j would require the computation, at each time step of the zeros of the approximate 
electrostatic field, which would increase substantially the cost, taking into account that we 
use high order approximations. Nevertheless, as we show here, the slight variation in the 
schemes does not affect the optimal accuracy of the methods. Furthermore, since the new 
definition of the fiux is still consistent, the mass and energy conservation can still be guar- 
anteed (even at the theoretical level). Also, here we demonstrate numerically that for the 
energy preserving scheme given in [4j , it is indeed necessary (and not a technical restriction 
due to the proof) to use finite element spaces spaces containing all quadratic polynomials. 

For the time discretization we stick to a simple fourth order explicit Runge Kutta (RK) 
method, the so-called RK4 or classic Runge-Kutta [39j. The reason for not using total 
variation diminishing (TVD) RK integrator is twofold. On the one hand, in our simulations 
we have observed no numerical evidence of any essential benefit of the TVD integrator over 
the standard RK method (probably due to the smoothness of the solution). On the other 
hand, since we focus on high order methods (for the space discretization), the time integration 
should be accomplished also with some high order time integration scheme. As is well known 
[37j, a fourth (or higher) order TVD RK, would require for the computation of the internal 
stages, the evaluation of the operator and its adjoint, due to the presence of some negative 
coefficients in the corresponding TVD-RK tableau. This would substantially increase the 
cost (and storage) of the overall procedure, without any significant benefit. With the fully 
discrete schemes, we verify numerically the theory developed in [4j; both the error analysis 
together with the conservation properties. 

The second goal of the paper is to validate the methods by studying their performance 
in approximating some of the classical benchmark problems in plasma physics. Here we 
consider the linear and nonlinear Landau-damping together with two benchmarks related to 
the two stream instability problems. We compare our numerical results with those available 
in literature, getting always at least the same outcomes. In particular, we show the benefit 
of using the energy preserving high order DG method for the numerical simulations (since 
no extremely refined meshes are needed and the code can be parallelized). 
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In the last part of the paper, we consider the apphcation of the schemes for the boundary 
value problem of the VP system studied in [30]. This problem models the evolution of a 
collisionless electron gas under the influence of a electrostatic fleld E in an interval [0, 1], 
with electrons emitted at one end and absorbed at the other end of the interval. Due to 
the absorbing boundary condition, it has been proved theoretically the distribution function 
/ might become discontinuous in flnite time, depending on the sign and magnitude of the 
electrostatic fleld at the boundary. Although the DG methods we consider in this paper 
were not originally designed to approximate such problem, we study here the ability of the 
methods to capture the discontinuity. The results however, are not completely satisfactory, 
since we do not always (at all the times) capture the behaviour of the solution predicted 
by the theory developed in [30j. A possible reason is the weak nature of the singularity, 
but it might also happen that as the time evolves the full discretization is adding too much 
artiflcial viscosity, which does not allow the methods to capture completely the singularity. 
This issue together with the tuning of the schemes to capture correctly the singularities (at 
all times) will be the subject of future study. 

The outline of the paper is as follows. In Section [T] we describe the main properties 
of the continuous problem and introduce the basic notations related to the discrete DG 
methods. We then introduce the numerical methods we consider discussing also their main 
properties in Section [2} In Section [3] we consider the full discretization and deal with the 
implementation issues related to the schemes. Section [4] is devoted to the validation and 
convergence study of the schemes. We present extensive numerical tests and consider the 
application of the methods for the simulation of Landau damping and two different tests 
related to the nonlinear two stream instability. In section [5] we examine the application of 
the considered DG methods for approximating a Vlasov-Poisson boundary value problem 
(no periodic boundary conditions). Finally, we derive some conclusion in section [6] 

Notation: Throughout this paper, we use the standard notation for Sobolev spaces (see 
[Ij). For a bounded domain 5 C M^, we denote by H^{B) the standard Sobolev space of 
order m > and by || • \\m,B and | • \m^B the usual Sobolev norm and seminorm, respectively. 
For m = 0, we write L'^{B) instead of H^{B). We shall denote by H^{I)/W the quotient 
space consisting of equivalence classes of elements of H^(X) differing by constants; for m = 
it is denoted by L^(X)/]R. We shall indicate by Ll{B) the space of L'^{B) functions having 
zero average over B. 

1. The Vlasov-Poisson system and basic notation 

In this section we introduce the Vlasov Poisson system and recall some of its properties. 
In the last part of the section, we also introduce the basic notation required for describing 
the numerical methods we consider. 

1.1. Continuous problem: the Vlasov-Poisson system. We consider a noncollisional 
plasma of charged particles (electrons and ions). For simplicity, we assume that the properties 
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of the plasma are one dimensional and we take into account only the electrostatic forces, thus 
neglecting the electromagnetic effects. We denote by / = f{x^v^t) the electron distribution 
function and by E{x^t) = the electrostatic field. The Vlasov-Poisson equations of 

the plasma in dimensionless variables can be rewritten as, 

df df df 

(1.1) + = ° ^ f^. x m x [o,t;], 

(1.2) = (^^t) e n,x[o,tfi 

where v denotes the velocity of the charged particles and p(x, t) is the charge density defined 
by 

p{x,t) = / f{x,v,t)dv V (x,t) X [0,t/]. 
Jr 

Let /o denote a given initial distribution /(x, 0) = /o(^, v) in (x, 'u) G [0, 1] x M. We impose 
periodic boundary conditions on x for the transport equation ( |l.lD , 

/(0,t-,t) = /(l,t-,t) V(t;,t)GEx[0,t/], 

and also for the Poisson equation ( |1.2D ; i.e., 

(1.3) $(o,t) = $(i,t), Vte[o,t/]. 

To ensure the well-posedness of the Poisson problem we add the compatibility (or normaliz- 
ing) condition 

(1.4) / p{x,t)dx= [ [ f{x,v,t)dvdx = l, VtG[0,t/], 
Jo Jo Jr 

which is the condition for total charge neutrality. To guarantee the uniqueness of its solution 
$ (otherwise is determined only up to a constant), we fix the value of $ at a point. We set 

(1.5) $(0,t) = VtG[0,t/]. 



Notice that (1.4) express that the total charge of the system is preserved in time. 

Through the paper we are only concerned with compactly supported solutions / of problem 



(1.6)-(1.2). We assume that a bounded set C M such that 

supp(/(x,^,0)) Usupp(/(x,^,t)) cn^xQy , Vt G [0,t/] , 
and so the Vlasov equation ( |1.1 ), can be (and will be) regarded in x [0,^/]: 
df df df 

(1.6) ^ + ^^'""'^^ ^ nxxn,x[o,tf]. 

The charge density is accordingly defined by 



(1.7) p{x,t) = f{x,v,t)dv V (x,t) G X [0,t/]. 
We define the total energy of the system as 

(1.8) £{t)= I f{x,v,t)^^dxdv+ I ]^\^x{x,t)\^dx VtG[0,t/]. 

The first term in the above definition represents the kinetic energy; the second, the potential 
energy of the system. 
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1.1.1. Properties. The Vlasov-Poisson system preserves in time many physical observables. 
We now briefly revise some: 



• Mass conservation: as already mentioned, the total charge of the system is preserved: 

(1.9) — / f{x,v,t)dxdv = VtG[0,t/]. 

dt Jq 

• LP -conservation: noting that div^^^^ (['U, — = one can deduce straightaway 
the conservation of all L^-norms of the distribution function: 

(1.10) / \\f{x,v,t)\\Pdxdv = VtG[0,t/]. 

dt Jq 

We will be particularly concerned with p = 1,2. 

• Total Energy: Following one can also show the following energy a-priori estimate: 



(1.11) ^'^W = ^(^/(^'^'^)^^^^^ + ^ ^\Eix,t)\^dx^ =0 yte[o,tf], 

where we have already used the deflnition of the electrostatic fleld E{x^t) = ^x{x^t) 
(compare with (1.8)). 



1.2. Basic notation and preliminaries for the numerical methods. Let {Th} be a 
family of partitions of our computational/physical domain Q = x = Qj. x [— L, L], 
which we assume to be regular [20j and made of rectangles. Each cartesian mesh Th is deflned 
as 

where 

h = [^i-i/2,^i+i/2] Vi = 1, . . . , A/"^; Jj = [^j-i/2,^j+i/2] Vj = 1, . . . , A/"^ . 
The mesh sizes hx and hy relative to the partition are deflned as 

< hx= max := - Xi-1/2, < hy = max /i^ := - Vi_i/2 , 

with and denoting the cell lengths of li and Jj, respectively. The mesh size of the 
partition is deflned as /i = max (/i^^, /i^). The shape regularity assumption implies that 
3ci, C2 > constants independent of h such that Cihy < hx < C2hy. 

We assume that = corresponds to a node of the partition along the '^-axis, i.e., Vj-1/2 = 
for some j in the partition of Qy = [— L, L]. We denote by {Th} the family of partitions of 
the interval Qx' ^h '= { li ' 1 ^ ^ ^ }• 

For /c > 1, let P^(/i) be the space of polynomials of degree up to fc, and let Q^(T^j) be 
the space of polynomials of degree at most k in each variable {{x^v)). We deflne the flnite 
element spaces: 

(1.12) V,' = {ijeL\l) : ijeF^h), yi,,i = i,...N,,}, 

(1.13) Zj; := {^eL'in): ^ e Q'=(r,,), VT,, = /, x J,, V^ ,i } , 

(1.14) W}: = {xeC°(X) : xeP'^(/i), V/i,^ = l,...iV,,}nL2(X)/M. 
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As is usual in the DG methods, we now introduce the the trace operators. We denote by 
(^h)t+i/2,v (^^)i+i/2,^ values of Lph at (2:^+1/2, v) from the right ceU x Jj and from 
the left cell li x Jj, respectively; 



for all {x^v) X J or in short-hand notation 

(1-15) (^/^)^tl/2,. = M^%i/2^ ^) ' (^/^)x,j+l/2 = ^h{x, Vf_^y^) , 



for all (x, v) G liX Jj. The jump [ • ] and average {•} trace operators of at '^), V'^ G 
Jj are defined by 

l^h\^+l/2^v := {^h)t+l/2,v - {^h)^+l/2,v ^ , 

Wh}i+l/2,v '= ^ i^h)t+l/2,v + (^/^)i"+l/2,^ ^ • 

For > 0, let : L^{X) — > Vj^ be the standard L^- orthogonal projection onto the finite 
element space V/^ defined locally, i.e., for each I < i < 

(1.16) / {P'{w)-w)qhdx = Vg^GP'(/,). 
By definition the projection is stable in L'^{I) 

(1.17) WP'Hhm,) < Mmi) e l\x). 

We denote by Vh ' L^{^) — > Z^ the corresponding two dimensional L^-orthogonal projec- 
tion; defined by Vh{w) = {P^ ® P^){w)] i.e., for all i and 

(1.18) / / {Vh{w{x,v))-w{x,v))iph{x,v)dvdx = {) V(^/,gP'(/,)®P'(J,). 
Also from its definition, its L^-stability follows immediately. 



2. Discontinuous Galerkin methods for the Vlasov-Poisson system: 

semi-discrete methods 

In this section, we introduce the DG methods we consider for approximating the Vlasov- 
Poisson system. The methods are those proposed in [U [5], but with some slight variation 
required for practical computations. Following f4^^ ^ we first describe the schemes for the 
Vlasov equation, proposing several options to modify the methods in [4j so that they allow 
for a feasible implementation. We then discuss the approximation of the Poisson problem 
(again following closely j^)- We close the section by discussing the main properties of the 
introduced schemes. Throughout the whole section we focus on the space discretization. 
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2.1. Discontinuous Galerkin approximation to the Vlasov equation. We now des- 
cribe the DG methods we consider to approximate the Vlasov equation ( |1.6D . For the time 
being, we assume we are given a finite element (conforming or nonconforming) approxima- 
tion of degree r to the electrostatic field E{x^t) = t), which we denote by G Wh- By 
El we refer to its restriction to The properties and characterization of E^ are discussed 
in next subsection. 

We denote by fh{0) = Vh{f{x^v^O)) the approximation to the initial data /(x,'u,0) com- 
puted using the orthogonal L^-projection onto the space Z^. Since the Vlasov equation is 
a transport equation, we construct the DG method in the usual way: given E^ G Wh find 
fh • [0,t/] — y such that 

(2.1) 5^5]4.(£;,;A,^,) = V^,GZ^ 

i=l 3=1 

where the bilinear form B^j{Eh] //j, Lph) is defined for each i, j and (fh G Zj^ as: 

IT dt 

(2.2) +/ ((vfh)(fiu)i+^/9,^ - ((vh)v>t)i-^/9.^ dv 



Bij{Eh;fh,^h) = f ^(phdvdx- [ vfh^^dvdx+ [ E^fJ^dvdx 



T dx It 



dv 



In (2.2) we have used the short hand notation given in ( |1.15 ). The numerical fiuxes are 
defined using the upwind fiux: 



(2.3) vfh 



V if w < 



vfh = {vfh} 



[fh], 



(2.4) Elfh 



Elf+ iiV'{El)>0 
Elfi; iiV'{El)<Q 
At the boundary dVt^ the numerical fiuxes are taken as 



El 



Elh = {EM + sign {r'{El)) . ^[A: 



{vfh)l/2,v = {vfh)N^+l/2,v^ {Eifh)x,l/2 = {Eifh)x,N.+l/2 = 0, V (x, ^) G X X J, 

SO that the periodicity in x and the compactness in v are refiected. Note that the numerical 
fiuxes as defined in ( 2. 3D and (2.4) are consistent. 



Observe that, unlike in [4, 5j the definition (2.4) of the upwind fiux E^fh involves a condition 
on the sign(P°(£'^)), rather than on sign(£'^). Obviously if Ej^ does not vanish inside li (and 
so it does not change sign inside /^), sign(£'^) = sign {V^{El)) and the classical definition of 
the upwind fiux is recovered. However, since E^ is a piecewise polynomial of degree k + 1 
approximation to the electrostatic field, it will in general change sign in some elements li of 
the partition I^- The classical definition would require to construct (at each time step) a 
partition of that adapts to the changes of sign of E^ by locating the zeros of E^ at nodes 
of the desired partition. Such process, although feasible in one dimension, might become too 
expensive and complicate unnecessarily the whole solution method for the Valsov-Poisson 



system (specially in higher dimensions). The definition (2.4) is considered for computational 
purposes. It allows to avoid computing the zeros of and re-meshing, at each time step, 
the partition 1^. 

For our computations of the one-dimensional problem we have also examined two other 



variants of the numerical fiux E^fh defined in (2.4), that do not require re-meshing, although 
they require a control on the sign of E^. We also show how this control on the sign(£'/j) can 
be done efficiently. The first variant we consider is given by 



iiEl>0 



(2.5) 



E'hfh 



Ehfh 

{EiU}^ [i3x*eli 
where we have used the weighed average 
(2.6) {El = ^^Elf+ + uj-Eir, 



UJ 



such that ^^(x*) = 

, UJ~ G [0,1] = 1. 



The parameter uj should be chosen so that the amount of upwind is tuned. Although based 
on heuristics, in our computations we have found that a good choice is given by 



(2.7) 



UJ 





max/. 


Eh 




max/. 


Eh 


+ 


min/, Eh\ 



UJ 





min/, Eh 




max/^ 


Eh 


+ 


min/, Eh\ 



The definition of the numerical fiux ( |2.5[ ) can be rewritten in the compact form: 

if ^ X* G li such that 
if 3 X* G Ij such that 



(2.8) Elfh 



{Elh] + ^ 
1, 



{Elh} + El{uj^--)lfh, 



Elix*) 
Elix*) 








Observe that this definition of the numerical flux is also consistent. 
The last variant we consider is defined by: 



(2.9) Elfh = { 



El fit 
Elfh 

r%Ei)f+ 

V'iEDf^ 



iiEi>Q 
ifEi<0 



if V^iEi) > and 3 X* e h such that EHx*) = 
if V^iEl) < and 3 a;* G li such that EHx*) = 0. 



Note that the numerical flux defined above in (2.9) is not consistent (it fails to be consistent 



in those elements where El changes sign, where we commit an error of order Oih)). 



Notice that both the weighted average approach (2.5) and the last definition ( |2.9[ ) require 
the knowledge of those elements of the partition Xh where Eh vanishes. This information can 
be obtained very easily (at least in one dimension), by checking the sign of coefficients of Eh 
expanded in a basis with Bernstein polynomial^ See the Appendix |A] for further details on 
Bernstein polynomials and how the detection of change of sign is implemented. 



^Bernstein polynomials are non-negative at everypoint of their domain 
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Mimicking ( 1.7[ ), we define the discrete density, ph{x,t): 



(2.10) ph{x,t) 



fh{x,v,t)dv 



J 



J2 I fhix,v,t)dv yx ex, Wt e [o,tf]. 



One of the nice properties of DG schemes, is that the conservation of the total mass is 
satisfied by construction. Next Lemma guarantees that the DG scheme (2.1)-(2.2) with 
fiuxes (2.3) and either (2.4) or (2.5) preserves the total charge: 

Lemma 2.1. Particle or Mass Conservation: Let fh G C^{[0,tf]; Z^), with k > 0, be 
the DG approximation to f , satisfying (2.1)-(2.2)^ with numerical fluxes defined as in (2.3) 
and either (|Z4|) or or (|2^. Then, for allt G [0,t/], 

(2.11) ^/ fh{t)dvdx = J2 fh{0)dvdx = J2 fodvdx = l. 



Although standard, we provide here the proof of the above Lemma for the sake of com- 
pleteness. 



Proof. {Proof of Lemma 2.1). Since fh{^) = ^h{fo) it follows from the mass conservation of 
the continuous VP system ( |1.4 ) and the definition of the L^-projection (1.18) that 



(2.12) J2 [ fh{0)dvdx = J2 f Vh{fo)dvdx = J2 f fodvdx = l. 



Now, let Tij be any arbitrary but fixed element in Th- By setting in (2.2) (/^/^ = 1 in Tij and 
(fh = elsewhere we find. 



Bij{EhJhA)=f^ I fhdvdx+ I [{vfh 



[(^lA).,,+i/2 - (^lA).,,-i/2] dx , 



dv 



where we have already used that such obviously satisfies {(ph) 



i+l/2,v 



i-l/2,v 



1 at 



the boundaries of Tij. Now, the above equation obviously holds for any i, since the choice 
of Tij was arbitrary. Therefore summing over all i and j the above equation, the fiux terms 
telescope and there is no boundary term left because of the periodic (for i) and compactly 
supported (for j) boundary conditions. Substitution now in (2.1) gives 

= J]Si,(E,;A,l) = - / fhdvdx = 0, 



and so integrating in time and using (2.12) we reach (2.11). 



□ 



2.2. Finite element approximation of the electrostatic field E, We now describe the 
methods we consider for approximating the electrostatic field E{x^t) = ^x{x^t). The discrete 
Poisson problem reads: 

(2.13) {^h)xx = l-Ph xeQx, $/,(l,t) = $/,(0,t). 
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The well posedness of the above discrete problem is guaranteed by (2.11) from Lemma 2.1 
which in particular implies 



(2.14) 



($/,),(l,t) = ($/,),(0,t). 



To ensure the uniqueness of the solution we also set $/i(0, t) = 0. 

Since in the Vlasov equation, the transport depends on to approximate (2.13) we consider 
a mixed finite element approach. For that purpose, we first rewrite problem (2.13) as a first 
order system: 



(2.15) 



Eh 



dx 



dx 



ph-1 X 



with boundary condition $/i(0,t) = $/^(l,t) = 0. We consider the following methods: 



2.2.1. Mixed Finite element approximation: we consider the one-dimensional version of Raviart- 
Thomas elements, KT^ k >1 [441 [T4] . In ID the mixed finite element spaces turn out to be 
the {W^^^,Vj^)-finite element spaces. Note that in particular, -^{W^^^) = V/^. For k > 
the scheme reads: find {Eh, ^h) ^ ^h^^ x ^^^^ ^^^^ 

(2.16) J Ehzdx + J ^hZa:dx = \/zeW^^\ 

(2.17) - J^{Eh),pdx = j^{ph - l)pdx yp G V^. 

We also refer to [6j , where the lowest order case was studied for the one-dimensional Poisson 
problem. 



2.2.2. Local Discontinuous Galerkin (LDG) method: The DG approximation to the first or- 
der system ( |2.15[ ) reads: find {E^, $/,) G V^^^ x Vj^^^ such that for aU i: 



(2.18) 
(2.19) 



EhZdx = - ^hZx dx + [{^hZ )i+l/2 - {^hZ^)i-l/2] 



WzeV^^\ 



EhPx dx 



•1/2 



{EhP^)^-l/2 = {Ph- l)pdx yp G V/:^' . 



(2.20) 



The numerical fiuxes {^h)i-i/2 and {Eh)i-i/2 for the LDG method are defined by: 

(^)^-l/2 = {^h}z-l/2 - Cl2[*/^l^-l/2 , 

{Eh)i-l/2 = {Eh}i-l/2 + Ci2iEh]i-l/2 + Ciil^hh-l/2 , 

where Cn = c{k + l^h'^ and |ci2| = 1/2. At the boundary nodes due to periodicity in x we 
impose 

{^h)l/2 = ($/i)7V^+l/2, {Eh)l/2 = (^/z)iV,+l/2. 

The method was first introduced in [21] for a time dependent convection diffusion prob- 
lem with Cii = 0(1). For the Poisson problem it has been considered in [18j in the one- 
dimensional case, and in [17j for higher dimensions. 
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2.2.3. Energy preserving LDG method LDG(v): We consider the DG approximation as given 
in (2.18)- (2.19), with numerical fluxes deflned by: 



{{^h)i-l/2 = {^h}i-l/2 - I$/*li-l/2 , 
• (v) 
{Eh)i-i/2 = {Ehh-1/2 + ^-^^Y^lEhji.y2 + Cul^hh-1/2 , 

with cii chosen as before, i.e., Cu = c(/c + l)^/i~^. Note that the above method requires the 
solution of two Poisson problems; one for ^ > and one for v < 0. Being one dimensional, 
this can be efficiently done without increasing the overall cost of the computation. 
This choice of numerical fluxes was introduced in [4J and extended to multidimensions in 
[5]. In both works it was shown that when combined with the classical standard upwind DG 
approximation for the Vlasov equation, the resulting semi-discrete DG method conserves the 
discrete total energy of the system. 



2.3. Properties of the numerical methods. We now briefly comment on the properties 
of the DG methods presented. The methods presented contain a small modiflcation of the 
schemes introduced and analyzed in [H [5] in a few elements. More precisely, the numerical 



flux Ej^fh in the DG scheme for the Vlasov equation has been redeflned through (2.4) and 
( 2.5[ ) (to allow for practical computations), in a few elements. Therefore, we expect that the 



resulting schemes will show similar stability, conservation and approximation properties. 



• Mass conservation: As we already showed in Lemma [2T} the total charge of the system 
is conserved. 

• L^-stabihty: we now comment on the L^-stability of the methods. In "5] the authors 
prove L^-stability for the DG schemes proposed there. Here, due to the modiflcation of 

the numerical flux (Ej^fh) in those elements where E^ might change sign, one cannot prove 
L^-stability (or at least the usual proof will not go through). Still, since the flux is modifled 
only in a few elements, it is reasonable to expect the methods to behave as if they were 
L^-stable. This will be verifled in the numerical experiments section [i} 

• Energy conservation: Finally, we deflne the discrete total energy: 
(2.22) 

£h{t)= / ^-h{x,v,t)dvdx+ / -\Eh{x,t)\^dx+ ^ ^ ' ^[$fe]?+i/2 VtG[0,t/]. 



Next result shows that also for the methods considered here with the modifled fluxes (2.4) 



and (2.5), when they are combined with the LDG(v) method for approximating the Poisson 



problem, the total energy of the Vlasov Poisson system is preserved. 

Theorem 2.2 (Energy conservation). Let k > 2 and let {{Eh.^h), fh) be the LDG(v)-DG 
approximation belonging to C-^([0, T]; (V^^^ x V^^) x Z^) of the Vlasov-Poisson syste m 
(1.2)^ where fh is the solution of ( |2.lD ^ (2.2)^ with numerical fluxes ( |2.3D and either ( |2.4[ ) or 

12 



(2.5), and the approximation (^h,Eh) solves of (2.18)- (2.19), with numerical fluxes (2.21). 
Then, the total discrete energy is conserved in time, 



(2.23) 



E / \vf Mt)dvdz + J2 I En{tfdx + c,^Yj^hmti/2] = 0. 



Proof. The proof follows exactly the same steps as the proof of [4", Theorem 5.1]. All the 

arguments used there carried over for the methods given here with the numerical flux E^fh 
modified as in (2.4) and (2.5). The reason is that the definitions of the fiuxes (2.4) and (2.5) 

are consistent, which is the only property needed for E^fh to ensure the conservation of the 
total energy. We omit the details for the sake of conciseness. □ 



As already noticed in [4^i5j, Theorem 2.2 requires the use of polynomial degree k >2. We 
will show in the numerical experiments that this restriction is not technical, but it is indeed 
required in practice. 



3. Fully discrete method and implementation details 



In this section we describe the time integration we consider and we discuss the details on 
the final solution algorithm. In last part of the section we also comment on the implemen- 
tation of the algorithm. 



3.1. Time integration. The DG methods presented so far are semi-discrete. For the time 
discretization we consider a simple fourth order explicit Runge Kutta (RK) method, the 
so-called RK4 or classic Runge-Kutta ^39j. While for conservation laws and other general 
nonlinear hyperbolic problems, in order to have good resolution of the shocks and disconti- 
nuities a total variation diminishing (TVD) RK should be used, here we have observed no 
significant differences (see Fig. 14), probably due to the inherent smoothness of the solution. 
Moreover, since we are concerned with high order methods (in space), the time integration 
should be accomplished also with some high order time integration scheme. As is well known, 
a fourth order TVD RK, would require for the computation of the internal stages, the evalu- 
ation of the operator and its adjoint, due to the presence of some negative coefficients in the 
corresponding TVD-RK tableau. Therefore, the cost (and storage) of the overall procedure 
would substantially increase, and from the experiments carried out, we have no numerical 
evidence of any essential benefit. This issue deserves surely a further theoretical study, that 
we plan to do in the future. For these reasons, although in general, is much safer to use a 
TVD Runge-Kutta method for solving hyperbolic problems, we have stick to the classical 
fourth order RK for the simulation of the Vlasov-Poisson system. 



We now describe the actual implementation. For the Runge-Kutta integrator, we take a 
uniform partition of the time interval [0,t/], {0 = to < < • • • < < • • • ^A^t = tf} with 
time step At := t^+i — t^. The RK4 solver updates the current approximate solution, from 
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fh{tn) to fh{tn+i)^ in four internal stages. The method is fourth order accurate in time, so we 
expect the errors coming from the time discretization does not affect much the ones coming 
from the space discretization. After choosing basis functions in (2.1)-(2.2) we arrive to a 
system of ODE's. 



(3.1) 



M^ = £(f,,E,,t), 



where M denotes the mass matrix, and f/^ and E/^ are the vector representation of the 
unknowns fh and E^^ respectively, in the chosen basis. The vector of unknowns, f/^, is arranged 
so that those degrees of freedom corresponding to the same element, say T^j, are in the same 
block. Hence, the structure of f/^ looks as follows: 



The mass matrix, M, is then block diagonal (in Lagrange basis) with each block of size 
(/c + 1) • (fc + 1), corresponding to the degrees of freedom in each element Tij. Therefore, 



the application of the RK4 solver to |3.1| can be done elementwise since it is only needed the 
inversion of the local mass matrices of size {k + 1) • {k + 1) where A; is a moderate integer, 
i.e. k < 12. This property allows performing the time marching from to t^^^ in parallel 
as we describe later. The local matrix inversion is done before starting the time integration 
and it is stored and saved for re-use in the whole computation. 

To advance in time from (r,f^,E;j) to (^+^ f^+\ E;J+^) the RK4 method proceeds in 4- 
stages: 



ki ■.= AtM-^C{{j:,El,t, 



k2 := At M-i£(f- + ki/2, E;:+'/^ + At/2) 
k3 := At M-i£(f^- + k2/2, E;:+'/^ t^ + At/2) 
k4 := At M-i£(f^- + ks, E;:+'/^ t^ + At) 



f^"+i := fj: + |(ki + 2k2 + 2k3 + k4) 

where 

E^+'/^ := SolvePoissonUsing[p^+^/^] = 1, ...3, 

and is computed using the corresponding approximate solution at that stage. 

Since Vlasov equation is a 2-dimensional transport problem, we choose At such that 

At ocmin (/i^ , h]) / m^x{\\Eh\\L^(n^),2L} , 

In the actual implementation we have tuned the At to ensure that the errors coming from 
the time discretization do not pollute the errors from the space discretization (that we want 
to observe) . In the experiments we have also used the TVD RK2 time integrator [371 SB] , 
to compare the results obtained with the two time integrators, and to see if there is possible 
advantage. (See the experiments on the Nonlinear Landau-Damping in Section [i]). It is the 
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following second order RK scheme: 

ki ■.= q + AtM-^£{fJi,El,t„ 



(3.2) 



where 

(3.3) E;;+^/^ := SolvePoissonUsing[p^+^/2], 

with the same expression for At. We have written it in the above form to evidence that the 
time advancing is done by means of linear combination of intermediate Euler steps (see [37j 
for further details). 



3.2. Solution algorithm. We now provide the pseudo-algorithm used in our computations 
to solve for {{fl, Ej:)}n>o given /o = v, 0)). 



(1) Given t/, At, /(x,'^,0), k (degree of polynomial), A^^^ and A^^ 

(a) Build a 2D mesh suitable for DG scheme using N^^ Ny^k (Vlasov part) 

(b) Build a ID mesh suitable for DG scheme using A^^^, {k + 1) (Poisson part) 

(c) Initialize the = Vh{f{x, v, 0)) 

(2) at each time: t^, Vn = 0, A^t 

(a) at each stage of RK4: 5 = 1, 4 

(i) Solve the Poisson problem and computing 

(ii) compute using C (/^(^-^)/^ t,) 

(b) compute the fj^^^ using {kg}. 

(3) Give the approximate solution at tf : fh{tf) 



Since the Poisson problem is one dimensional, its solution (Step 2. a (i)) is done using an 
exact solver. This is pre-computed at the initial time step, and then at each Step 2. a (i), only 
two matrix vector multiplication are required. For higher dimensions, such solution process 
should be done iteratively and with an appropriate preconditioner [36l [2] . We should mention 
during actual implementation and verification of the code, we have observed that in order to 
guarantee accuracy and to ensure the conservation of mass and the total energy, it is essential 
in every time step to compute at each stage of the RK method. If on the contrary, one 
uses EJl at the all stages involved in the evolution from tn to tri+i, the high order accuracy 
and conservation properties of the methods are lost. 



3.3. Comments on the implementation. Now we describe a few issues related to the 
practical implementation of the scheme. As mentioned before, the fact of using an explicit 
ODE solver for (3.1) together with the block diagonal structure of the mass matrix M allows 
to perform time marching in a parallel manner. Here we briefiy demonstrate the structure 
of parallelization algorithm for time marching which is very useful when the size of the semi- 
discrete form becomes very big. Denoting the number of computational threads by A'cpu? we 
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partition the elements in Th into {Th,m} for m = 1, A^cpu where 



7" I lA^cpu 

Moreover in the same way we decompose the vector of degrees of freedom fh into {fh,m} and 
kg into {ks,m} for m = 1, A^cpu- Then for the paraUel time marching we have 



(1) Given A^cpu (number of threads), {Th,m} (partition of Th) 

(2) at each time: tn^^n = 0^ ...^ Nt 

(a) at each stage of RK4: 5 = 1, 4 

(i) Solve the Poisson problem and compute 

(ii) on each partition Th^rn' m = 1, A^cpu 

(A) compute fc,,^ using £ (/r^^"'^^', C^^^"'^^', 

(iii) construct kg = U^^^^kg^rn • 

(b) compute the fj^^^ using {kg}. 

(3) Give the approximate solution at tf : fh{tf) 



Note that Step 2.a.ii. is the parallel part of the algorithm. The structure of the above 
algorithm is well-suited for OpenMP library for parallelization. We have run the algorithm 
on maximum 8 thread using OpenMP when the size of system becomes very big. In order to 
boost the number of threads one should modify the algorithm structure to run under MPI 
library which was not needed during our numerical experiments. 



4. Numerical experiments 



In this section we present several numerical test to assess the performance of the introduced 
methods and to validate their properties. We start with a convergence study using a forced 
VP system for which the exact solution can be explicitly computed. Then, we present the 
results obtained for simulations of some of the classical benchmark test for the VP system 
with periodic boundary conditions. More precisely the tests we consider are: 



(1) Convergence study: forced Vlasov-Poisson system, 

(2) ID weak Landau damping, 

(3) Nonlinear (strong) Landau damping, 

(4) Two stream instability I, 

(5) Two stream instability II. 



To study and asses the ability of the proposed DG schemes to conserve the physical invariants 
of the continuous Vlasov-Poisson system, we have computed for the different tests, the time 
evolution of the deviations from their initial values of the quantities that are conserved by 



the continuous Vlasov-Poisson system (see Section LI). The value of these quantities at time 
t = can be always computed using the discrete initial data (which is the L^-projection of 
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the continuous initial data as defined in ( 1.18D ). We will study the following deviations (from 
conservation) : 



\£hit)-8hm 

4(0) 

/^^ ph{t, s)ds - fh{x, V, 0)dxdv 
J^fhix,v,0)dxdv 

||A(t)IUi-||A(o)IUi 
IIA(o)IU^ 

Most of the computations are carried out with the RK4 as described in section[3 and the 
energy preserving DG-LDG(v) method, with numerical fiuxes (2.3)-(2.5) and (2.7) for the 
Vlasov discretization and (2.21) for the LDG(v) discretization of Poisson. When using the 
other described methods (in space and time) we will explicitly say. 
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Figure 1. Forced VP: Convergence diagram for the L^-error in the approx- 
imation of the distribution function, /i-convergence diagram (left) 
and fc-convergence diagram (right). 



4.1. Test 1: convergence for a forced Vlasov-Poisson system. The aim of this test 
is to validate and assess the convergence properties of the proposed DG schemes. We set 
= [— TT, tt] and Qy = [—4, 4] and consider the following VP system in Q = Qx x Q^: 



(4.5) 



ft + vfx - E{x, t)fy = V^(x, V, t) 
d 



dx 



E{x^ t) = p(x, t) — \p-K 



(x, v,t) ^Vtx 

X G [— TT, tt], 
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E{x,t) = ^sin(2x - 27rt) 



where the right hand side ip{x^v^t) is chosen so that the exact solution (f^E) of (4.5) 
(complemented with periodic be in x and compact support in v) is given by 

/(x, V, t) = {2- cos(2x - 27Tt)} e-i(^^-^)' (x, ^, t) G x M+, 

The forcing term ilj{x^v^t) in ( |4.5[ ) is defined by 

^, t) = e-i(4^-i)' ({(47^ + 2)^ - (27r + a/^)} sin(2x - 27rt) 
+ Vtt (1/4 - ^) sin(4x - 47rt)) . 
Observe that the solution (/, E) is periodic in time with period 1: 

f{x,v,l) = f{x,v,Q) \i{x,v)en. 
Therefore, we have performed the the computations up to time tf = 1. In Figure [T] are given 



I 




Figure 2. Forced VP: Convergence diagram for the error in the approxima- 
tion of the electrostatic field E^. 

the convergence diagrams for error of the distribution function approximated by the 



energy preserving DG method (2.18)-(2.19)-(2.21). The /i- convergence diagram is depicted 
on the left figure for several polynomial degrees = 1,2,...6. In agreement with the 
theory in [4], optimal order of convergence A; + 1 is achieved when using the method with 
approximation degree polynomial k. On the right figure we have also represented the k- 
convergence diagram (error versus polynomial degree) varying the polynomial degree k = 
1, 2, 4, 8, 12 for different meshes 20x20, 40x40 and 80x80. Although not proved theoretically, 
the results in the figure seem to indicate an exponential rate of convergence in this case. 
The corresponding convergence diagram for the error in the electrostatic field is given in 
Figure [2] Also in this case the theory in [4J is verified. Note that as predicted in [4J the same 
convergence rates are obtained for the approximations of fh and E^^ even if the approximation 
to Eh is done using polynomials one degree higher. Observe though that the error in E^ is 
much smaller. We now study the effect (in accuracy) of using the different Poisson solvers 
in the DG method for the Vlasov Poisson system ( |4.5[ ). Together with the LDG{v) which 
gives the energy preserving scheme, we have also run the computations using the classical 
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Ij KillKjl 
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1/20 
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3.0134 X 10-2 




1/40 


6.4640 X 10-3 
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6.4623 X 10-3 


2.2213132 


6.4623 X 10-3 


2.2213157 


1/80 
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3.0922407 


7.5775 X 10-4 


3.0922409 


fc = 3 


RT3 


LDG(w) 


LDG 


h 


L? error 


order 


L? error 
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L?' error 


order 


1/20 


5.8300 X 10-3 




5.8295 X 10-3 




5.8295 X 10-3 




1/40 


3.6364 X 10-"^ 


4.0029199 


3.6361 X 10-"^ 


4.0029180 


3.6361 X 10-4 


4.0029180 


1/80 


2.2582 X lO-'^ 


4.0092237 


2.2580 X lO-'^ 


4.0092260 


2.2580 X 10-^ 


4.0092260 



Table 1. Forced VP: Convergence rates and errors — //i(t/)||o,r^ for 

k = 2 (top) k = 3 (bottom) using different Poisson solvers (with 
k + 1 polynomial spaces). 



LDG (with numerical fluxes deflned in (2.20) and flnite element spaces {Vj^^^ ^Vj^^^) and 
the mixed finite element approximation (2.16)-(2.17) with the RT^ (Raviart-Thomas) finite 
element spaces {W^^'^^ ^h)- Table [l] are given the L^-errors ||/ — fh\\o,Th convergence 
rates for different mesh sizes and polynomial degrees k = 2 and k = 3. As it can be 
observed, all the methods considered seem to yield approximations with the same accuracy 
and convergence properties. We now study how the discrete mass and discrete energy are 
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mesh 40 X 40, k = 6 
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mesh 40 X 40, k = 6 
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Figure 3. Forced VP: Time evolution of the deviation from conservation of 
mass and total discrete energy 



preserved in time. In Fig. 3, are depicted the time evolution of the relative error of the mass 



conservation (left), and energy conservation (right). The corresponding diagram for (4.4) is 
given in Figure [4} The results are obtained on a uniform mesh 40 x 40, using the DG energy 
preserving method with polynomial degree k = 6. Notice that the errors in the graphics 
are close to machine precision, which indicates the ability of the method to conserve the 
properties of the system. 
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Figure 4. Forced VP: Deviation from conservation of the L^-norm of fh- 
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Figure 5. Weak Landau Damping: Time evolution of the amphtude of the 
electrostatic field. 

4.2. ID weak Landau damping. Typically, in most works in literature concerned with 
the simulation of Landau damping (linear, weak or nonlinear) for the VP system, the compu- 
tational domains for the phase space Q = x are set to Qx = [0, 47r] and Qy = [—5, 5]. 
However in our computations, we have found that using such f^, after some t > 0, the 
approximate distribution function would not be of compact support in v. More precisely 
we found that ^ 10~^ for large time which is far from 0, specially for high 

order accurate approximations as those considered in this paper. Therefore, to ensure the 
compact support in v of fh we have set = [—10, 10] for our computations. In this case, 
//^(x, 'z;, ^ 10~^^ which can be obviously regarded as 0. We take as initial data 

1 _2l 

(4.6) /(x, V, 0) = —= {1 + a cos(K x)) e 2 x e Qx ^ ^ j 

where Qx = [0, 47r] and = [—10, 10]. In ( |4.6| ), a is the size of the perturbation and K refers 
to the basic mode of the electrostatic field. Here, we have set K = 0.5 and a = 0.01 so that 
the perturbation is small and therefore the linear theory can be used. We have computed the 
approximate solution {E^^ fh) over a mesh 60 x 60, using the energy preserving DG method 
with polynomials of degree /c = 4. 

In Figure [5] we plot (in a semi-log diagram) the time evolution of the L^-norm of the 
electrostatic field Eh{t). As it can be observed from the graphic, the amplitude of the 
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Figure 6. Weak Landau Damping: Time evolution of the relative error 



(deviation) of the total energy (4.1)(right) and norm (4.3) (right) 
for the weak Landau damping. 



electrostatic field decreases exponentially in time up to some recurrence time T^^ after which 
it oscillates, in agreement with Landau linear theory. We have fitted the line (in log-scale) 
c exp(— 7t) at the local maximums ||£'/^(t)||o. The obtained damping rate 7 of the oscillations 
is 7 = —0.153272, which is in good agreement with those results found in literature (compare 
to -0.1533 in [45j). 



k - 4, a - 0.01 




90 100 



Figure 7. Weak Landau damping: Mass conservation 



In Figure [6] we have depicted the time evolution of the deviation from conservation of 
the discrete total energy £h{t) and the L^-norm of the approximate distribution function //^. 
The corresponding diagram for the mass conservation is depicted in Fig. [7| Observe that 
the deviation from conservation for all these quantities is close to machine precision, which 
assess the good conservation properties of the scheme. 

In [29] the authors study the performance of eulerian solvers based on finite volumes and 
finite differences , and in particular they compare the estimated with the theoretical time 
predicted from the free streaming case given by the formula = We have run some 
computations to see whether we could find a similar relation for the DG methods and so 
possibly depending on the polynomial degree k. The results are given in Fig. |8) However, 
from these results it seems to us difficult to provide a closed formula or relation for the DG 
methods. 
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Figure 8. Weak Landau Damping 



4.3. Nonlinear (strong) Landau damping. Nonlinear Landau damping is regularly used 
to assess the performance and properties of Vlasov-Poisson solvers (see [271 ED ESI ESI HD 
[33J. 



We take as initial data f{x^v^O) the function in (4.6), but we now set a larger amplitude 
of the initial perturbation of the density a = 0.5 and take K = 0.5. The computational 
domain is taken as for the weak case; Q = x with = [0,47r] and Qy = [—10, 10]. 
In this case, //i(^, ^ 10~^^ for large t, and so the compact support of fh in v is 
guaranteed. 



For this test, the Landau Linear theory cannot be applied since now the non-linear effects 
become important. However we will compare with other results obtained numerically in 
literature. 



In Figure Fig. 9 we plot the time evolution of the (log of the) ||£'/j(t)||o, computed with 
the DG-DG(v) method over a mesh 100 x 160 using polynomials of degree k = 3. On the 
right diagram in the same figure, we have represented the corresponding time evolution of 
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Figure 9. Non-linear Landau Damping: Time evolution of the amplitude 
of the electrostatic field. Estimated coefficients Ci = 2.279673, 71 = 
-0.292285, C2 = 0.015228 and 72 = 0.086126. 



mesh 


7decay 


^decay 


7growth 


^growth 


50 X 80 


-0.292286 


2.279682 


0.085114 


0.015669 


100 X 160 


-0.292285 


2.279673 


0.086126 


0.015228 


150 X 240 


-0.292285 


2.279673 


0.086116 


0.015232 



Table 2. Non-linear Landau Damping: 



Estimated values for coefficients 



of the fitted functions cexp(— 7^); different mesh size. 



the deviation from conservation of the discrete total energy of the system ( |4.lD . Observe 
that the amplitude of the electrostatic field decreases exponentially initially t G [0, 10] and 
then increases exponentially (for t G [20,40]) and after that it oscillates periodically. We 
have also depicted in the same diagram, the lines (in log-scale) cexp(— 7^) fitted at the local 
maximums of ||£'(t, .)||i^2 for t G [0,10] (initial decay) and for t G [20,40] (growth). The 
estimated coefficients and damping rates are Ci = 2.383814, 71 = —0.305920, C2 = 0.015360 
and 72 = 0.085241. They are in good agreement with numerical simulations presented in 
the literature: the estimated 7decay is same as the one reported in [45j —0.292; and close to 
the one obtained in [19]: —0.281. The growth rates are also in good agreement with those 
reported in literature (cf. [45]). We have also run this test using different meshes to study 
the effect of the mesh refinement and of increasing the polynomial degree k on the estimated 
values for the damping and increment rates. The estimated coefficients of the fitted functions 
cexp(— 7/;) are given in Table [2] (for the mesh refinement) and in Table [3] (for the increase in 
the polynomial degree). As expected, by refining the mesh or by increasing the polynomial 
degree, the estimated coefficients have less error. 



In Figure [T0| we have represented the corresponding time evolution of the deviation from 
conservation of the discrete total energy of the system (4.1). Note that up to t ^ 10 the 
deviation from conservation of total energy (quantity (4.1)) decreases to 10~^^ and after 
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mesh 100 X 160 



mesh 


Tdecay 


Cdecay 


Tgrowth 


^growth 


k = l 


-0.280475 


2.177524 


0.087209 


0.014650 


k = 2 


-0.292584 


2.244969 


0.085946 


0.015390 


k = 3 


-0.292285 


2.279673 


0.086126 


0.015228 



Table 3. Non-linear Landau Damping: Estimated values for coefficients 
of the fitted functions cexp(— 7t); different k. 



mesh 100 x 160, fc - 3 




t t 

Figure 10. Non-linear Landau Damping: Time evolution of the deviation 
from conservation of the total energy of the system. 



that due to the process of filamentation there is a slowly increment until t ^ AO when strong 
oscillations occur in ^'-direction. Therefore, the discrete total energy of the system is con- 
served with a relative error of order 10~^^, which to our knowledge has not been obtained 
before in literature. This property indicates that our scheme gives an accurate description of 
macroscopic values (physical quantities defined by the moments of the distribution function 
with respect to v). 



The evolution in time of (|4.3|) and (4.4) is depicted in Fig. 11 
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Figure 11. Non-linear Landau Damping: The evolution of 

1 1 //ill ^2 in a semi-log scale where using a mesh 100 x 160. 



In Figure [12] are represented in phase space the approximate distribution function fh at 
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(c) t = 45 
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Figure 12. Non-linear Landau Damping: approximate distribution func- 
tion at different times, computed over a mesh 100 x 160 with the 
DG-LDG(v) method with k = 3. 



different times, obtained with the DG-DG(v) method over a mesh 100 x 160 using polyno- 
mials of degree k = 3. From the figures it can be observed the detailed structure of the 
solution captured by the energy preserving DG method. To asses the ability of the scheme 



to capture the possible strong oscillations in 'u-direction, we have also plotted, in Figure [13 
the corresponding profiles in v of each fh{t)^ defined by: 

fh{x,v,t)dx . 



We have also performed the simulations with the RK2-TVD time integrator, to study how 
the use of such integrator affects the accuracy and conservation properties of the DG-DG(v) 
method. In Figure [Ml we have depicted the time evolution of the amplitude of the elec- 



trostatic field (left diagram) and the time evolution of (4.1) (right diagram), when using 
k = 3 over a mesh 100 x 160. These graphics should be compared with Figures [9] and [10 
(left), respectively. As regards the energy conservation, notice that the errors are higher 
than those obtained with the RK4, although this can be an effect of the lower order of the 
time integration used. In any case, the use of TVD integration does not seem to show any 
special advantage for this problem. 
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Figure 13. Non-linear Landau Damping: Profiles in v of the approxi- 
mate distribution function //^(^) at different times. 



mesh 100 x 160, k - 3, using RK2TVD 




mesh 100 x 160, k - 3, using RK2TVD 




Figure 14. Non-linear Landau Damping: results with RK2-TVD time 
integrator 



We now study and compare the conservation properties of the DG method when using the 
two approaches ( 2.4[ ) and (2.5) for defining the numerical fiux E^fh. In Figure |l5| are depicted 
the time evolution of the error in the total energy, obtained with both approaches for A; = 3 
(upper diagrams). The result obtained with the fiux defined as in (2.4) are represented on 
the left diagram. The one obtained with the weighted average modification given in (2.5) are 
given on the right diagrams. On the bottom diagram of the same figure is given the result 
corresponding to the use of the non-consistent definition (2.9). As can be appreciated from 
the results, the use of either (2.4) or (2.5) leads to the conservation of the total energy, while 
the definition (2.9) does not. These results confirm Theorem 2.2 Although not reported 
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here, the same set of experiments was run with k = 1,2. For the lowest order A; = 1, similar 
results were obtained for all the approaches, with no significant diflFerences among them. 



to ^ 
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mesh 100 x 160, k - 3, using signCP^iEl)) 
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Figure 15. Non-linear Landau damping: time evolution of the deviation 



(4.1) of the discrete total energy computed over a mesh 100 x 160 



with k = 3 with the different implementations of the numerical 

flux Ej^fh. Definition ( |2.4D (left upper diagram); definition ( |2.5D 
(right upper diagram), non-consistent definition ( |2.9[ ) (bottom di- 
agram) 



4.4. Two stream instability. This is a standard benchmark for checking the reliability of 
the schemes to face the strong oscillations. We have set Qx = [0, 47r] and Qy = [—10, 10] (as 
discussed in |4.3[ to ensure the approximation is compactly supported in Qy at all times of 
the computation) and take as initial data for the VP system ( 1.6[ )-( [L7 ), 

.,2 o 



(4.7) 



/(x,^,0) 



{2 - cos(K(x - 2^))} (x, v) e ^ 



with K = 0.5. In Figure 16 are represented the approximate solutions (in phase space) 
obtained at different times t = 15, 30, 45 and 60 with the energy preserving DG method 
using polynomial degree k = 3 over a mesh 150 x 150. The evolution in time of the deviation 



from conservation of the discrete total energy (4.1) and the L -norm of the solution (4.3) are 



depicted respectively, on the left and right diagrams in Figure 17, Notice that also in this 



case, the total energy of the system is preserved up to machine precision. 



We now compare the effect of using different Poisson solvers in the conservation proper- 
ties of the final DG scheme for the Vlasov Poisson system. In Figure 18 are depicted the 
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Figure 16. Two stream instability: solution of the VP system for two 
stream instability at different times using mesh 100 x 160 and 
k = 3. 



10-10 

10-11 
10-12 
D^- 10-13 

^ 10-1^ 
10-1'^ 



mesh 150 x 150, k = 3 




10 20 30 40 50 60 70 80 90 IOC 

t 



90 100 



Figure 17. Two stream instability: evolution of the relative error in 
and total energy in a semi-log scale 



time evolution of the deviation from conservation of total energy (4.1) obtained when using 
the mixed finite element approximation (2.16)-(2.17) with the RT^ (Raviart-Thomas) finite 
element spaces {W^^^ ^Vj^) (— — —); the classical LDG (with numerical fiuxes defined in 
(2.20) and finite element spaces {V^+\ V^+^) (- + -+) and the LDG(^;) method with fluxes 
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( 2. 3D (— o o ). On the left graphic are given the results obtained over a fixed mesh 

100 X 160 with = 1; on the right graphic those obtained with k = 3 and over a mesh 
40 X 40. From these graphics it can be appreciated that by using k = 1 the total energy is 
not conserved and the error committed is of 5 orders of magnitude higher than for A; = 3 
even for a fine mesh. Notice also that for A; = 3 even when using a very coarse mesh (40 x 40) 
the total discrete energy is conserved with error of order 0(10"''). 



mesh 100 x IQO, k -- 



I 




7 X 10-^ 



mesh 40 x 40, k = 3 



4 X 10-^ 




Figure 18. Two stream instability: time evolution of the deviation ( |4.lD 
for different Poisson solvers; the red line is obtained using LDG('^), 
the cross symbols are RT^ and plus symbols are classical LDG. 
k = 1 (left) and k = 3 (right). 



We now study the effect of mesh refinement for the conservation properties of the scheme. 
Here, we check the effect of refininging x and v separately. In Figure [19] are depicted the 
time evolution of (4.1) (top graphics), (4.3) (center) and (4.4) (bottom graphics) for k = 3. 
The graphics on the left show the effect of refinement of the x variable. On the right are 
given those obtained by refining only the v variable. From the graphics it can be observed 
that the and norms of fh depend more on refinement of the v variable but they seem 
to be insensitive to the refinement on the x direction. For the total discrete energy, however, 
the opposite effect can be appreciated from the top graphics. It seems to depend more on 
the refinement in the x variable (which is expected since the potential energy comes from 
the Poisson coupling). 



4.5. Two stream instability II. We set now Qx = [0, 137r] and Qy = [—8, 8] and take the 

initial data as in [43j 

(4.8) 

(1 + 0.05cos(Xx)) / / (^-'^)^\ / {v + wf 



/(x, V, 0) = = exp — + exp — {x,v) eft , 

with Vth = 0.3, w = 0.99 and K = j^. The initial data for this test consists of the two 
unstable flow moving in the opposite direction of each other. The approximate distribution 
functions fh obtained at time t = 70, with the DG-LDG(v) using a mesh 256 x 100 and 



A; = 1, 2, 3 are represented in Figure 20 



In Figure 21 are depicted the time evolution of the deviation from conservation of the total 
discrete energy. The left diagram shows the effect of mesh refinement in the quantity (4.1) 
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Figure 19. Two stream instability: time evolution of the deviation from 
conservation of the total energy (4.1) (top graphics), L^-norm of 



fh (4.3) (center) and norm of fh (bottom graphic). Diagrams 
on the left show refinement in x; on the right, refinement in v. 



80 90 




80 90 



80 90 



for the DG-LDG(v) method with k = 3. On the diagram on the right, we have represented 
the time evolution of such deviation for different polynomial degrees = 1, 2, 3 using a fixed 
mesh 256 x 100. From the right diagram it can be appreciated that k = 1 does not yield 
to the conservation of the total energy; in fact the error is almost 10 orders of magnitude 
higher than for > 2, even on a very fine mesh. This result indicates and confirms that the 
hypothesis A; > 2 in Theorem 2.2 is indeed necessary for the DG-LDG(v) scheme to preserve 
of the total energy of the system. 



5. Boundary value problem: a nonlinear plane diode 



In this section we consider the extension of the proposed DG schemes to approximate a 
Vlasov-Poisson boundary value problem. 
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Figure 20. Two stream instability II: solution at t = 70 using mesh 256 x 100. 
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Figure 21. Two stream instability II: time evolution of the deviation from 
its initial value of the total discrete energy. 



5.1. Model Problem: a nonlinear plane diode. We first describe the model problem, 
revising also some key results from [30j related to its theoretical analysis. Denoting by 
E '= the Vlasov-Poisson problem we consider is written as: 



(5.1) 
(5.2) 
(5.3) 



dxE = p 



fdv 



$(0,t) = 0, $(l,t) = A(t) 



ix,v,t) e [0,1] X [-L,L] X [0,tf], 
ix,t) e [0,1] X [0,tf], 
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with boundary conditions 



(5.4) 



f{0,v,t)=g{v,t), v>0 Vte[0,t/], 
f{l,v,t)=0, v<0 yte[o,tf], 



where g{v^t) is some given function. The system is complemented with an initial data 
/(x,'^,0) = fo{x^v). Note that the Dirichlet boundary conditions for the Poisson problem 
imply 



(5.5) / E{x,t)dx = X{t), VtG[0,t/]. 

Jo 

This problem models the evolution of a collisionless electron gas under the influence of the 
self-consistent electrostatic fleld E in the interval [0, 1], where electrons are emitted at one 
end of the interval {x = 0) and absorbed at the other end {x = 1). Due to the absorbing 
boundary condition, in general the distribution function / might become discontinuous in 
flnite time. Existence and uniqueness of solutions of nonlinear boundary value problems of 
this kind have been studied to a certain extent. The stationary one dimensional problem 
was considered in [31j and the higher dimensional case in j38l [29] . For the time dependent 
problem (5.1)-( 5^ -( |5.3 )-(5.4), weak solutions were constructed in [H^^^J. 



In [30j the authors carry out an study of the regularity of the distribution function solu- 



tion of the nonlinear plane diode (5.1)-(5.2)-(5.3)-(5.4). In particular they show that under 



rather general assumption on the data /o and that the electron distribution is of bounded 
variation (BV) as a function of x and v. This guarantees the uniqueness of the nonlinear 
plane diode. Furthermore, the authors study the influence of (the magnitude of) the exter- 
nal voltage A on the regularity of the solution. They show that if initially A(0) is large the 
solution might become discontinuous, while if A(0) vanishes or if it is small and the inflow 
boundary g{t^v) is non- vanishing, the electron distribution is of class C^{Vt) for all time. 

Prior to recall a couple of results from [30j, that allow for quantifying what is mean by 
"small" or "large" external voltage, we introduce some further notation. We denote by 75 
the singular set, 

(5.6) 7^ = {x = 0, ^ = 0} U {x = 1, v = 0}, 
and by 7J" the incoming set at the boundary x = 0: 

(5.7) 7o+ = {(0, ^, t) I ^ > 0, < t < tf} . 
We also deflne the norm: 

(5.8) |||[/o,^]||| := ||/o||ci + \\g\\c^ + Ib^x/olloo + ||^5^;/o||oo + \\v~^dtg\\oo + \\v~^d^g\\oo • 

Next result corresponds to [30l Lemma 4.3] 
Lemma 5.1. ( ]30l Lemma 4.3]J Let fo^g G C^{Q) with compact support be such that 

(5.9) TV[fo]+ [ {{l + v)\gv\ + \gt\) + \\v'fo\\oo + \\v'g\\oo<oo forsomep>2, 
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where 7J" is the incoming set at the boundary x = defined in (5.7). Assume they satisfy 
the following compatibility conditions: 



(5.10) 



fo{0,v)=g{v,0) v>0, 

fo{l,v) = v<0, 

dtg{v, 0) + vdJoiO, v) + E{0, 0)dJo{0, v) = v > 0, 

vdJo{l,v) + E{l,0)dJo{l,v) = < 0, 

and that |||[/o,fl']||| < co^ where |||[-]||| is the norm defined in (5.8). Further 

(5.11) /o(^5 7^ /o(0, 0) for (x, v) ^ (0, 0) and x, small. 

Then, if the external voltage satisfies 



(5.12) 

then f{x^v^t) is not continuous. 



A(0) > / / {l-x)fodvdx 
Jo J-L 



In absence of external voltage, the solution can be shown to be of class C^. 

Theorem 5.2. ([30, Theorem 4. 4]. J Let A(0) = 0. Assume fo^g G C^{Q) with compact 
support be such that \\\[fo^g]\\\ < oc and assume they do satisfy ( |5.9[ ) and ( |5.10 ). Then 

fix, V, t) e C an X [0, tf]) \ 75) nc\nx [o, t/]). 



5.2. Numerical schemes. To approximate the boundary value problem, (5.1)-(5.2[)-(5.3)- 



(5.4), we need to modify slightly the definition of the DG schemes given in Section^ The 



change of sign of the self-consistent electrostatic field in the Vlasov equation, forces the 
following change in the definition of the method: find (//^, Eh) : [0, tf] — > x such that 



where 



i=i j=i 



BijiEh, fh, ^h) = j -^^h dvdx- j vfh^^ dvdx- j Elfh^^ dv dx 



(5.13) 



+ 



h L 



dv 



xJ-1/2 



with the numerical flux vfh defined as in (2.3) and the definition of E^f is changed to 
(5.14) Ejjh=\ EjX = {Elh}-sign{V'{El)).^lh}. 
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Or, if using the weighted average 
(5.15) EiX = 



{Elfh}+^-^lfh\ 



with uj , uj^ 



1 — u~ defined as in (|2.7|). 



if$x* e k such that Ei{x*) = 
if3x* E k such that Ei{x*) = 



At the boundary dO,^ we still reflect the compactness in v in the numerical flux El^fh] 

while for the x-boundary nodes we account for the inflow boundary data: 

{vfh)i/2,v = v{ff; - g{v, t)) > 0, vej. 



To approximate the Poisson problem, we consider the LDG(v) method deflned in (2.18)- 



(2.19)-(2.21 ). The only modiflcation required is on the right hand side of the Poisson problem. 



since now there is no neutralizing background. At boundary nodes, the numerical fluxes are 
defined to account for the boundary conditions (|5.3[): 



(5.16) 



($/^)l/2 = ($/^)7V.+l/2 = A(t) 

(Eh)l/2 = {^/^}l/2 + Cll[$/,]l/2 , 
(^/i)iV,+l/2 = {^/i}iV,+l/2 + Cnl^h - A(t)]iV,+i/2 , 

where Cn = c{k + l)^/i~^. For the time discretization we use the RK integrator described in 
Section [3l 



5.3. Numerical experiments. We now present some numerical simulations obtained with 
the DG-LDG(v) method for the nonlinear plane diode . The main goal of this section is to 
verify if the proposed methods are able to detect and capture the singularity of the solution, 
when the external voltage Aq = A(0) is large and satisfies the condition (5.12). Following 
j, we consider the following initial condition 

fo{x,v) no{x)^=v^ exp{-v^/2) 



no{x) 



and the infiow boundary data: 

g{v,t) = 



271 



271 

(l + 7x)(l- 4x^)4 xe [0,0.5] 
else 

'exp(-^V2) Vt G [0,t/]. 



We have taken L = 10 to set the domain in velocity, since both /o and g are of compact 
support in [— L, L] = [—10,10]. It can be checked that the data satisfies the smoothness 
conditions |||[/o,5']||| < oo and (|5.9D together with the compatibility conditions given in 



(5.10). We have computed the approximate solution to (5.1)-(5.2)-(5.3)-(5.4) for different 
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Figure 22. Plane Diode: A(0) = 0. Evolution of ||VA||l2 (left) and 
(right) for different mesh size. 



values of the external voltage A(0). The main goal is to verify that the presented DG methods 
are able to detec t th e smoothness of the approximate solution. That is for A(0) = 0, in 
view of Theorem IK2| f{x,v,t) G {{^ x [0,t/]) \ 7^) nC^{n x [0,t/]), and we expect the 



approximate solution to reflect such regularity. 
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Figure 23. Plane Diode: Evolution of ||V//j||l2 for different mesh size (left) 
and different polynomial degree (right). 



In Figure 22 we have depicted the time evolution of L^-norm of Vfh and L^-norm of fh 
for different mesh sizes. As can be observed from the flgures, both quantities are flnite and 
bounded and converge toward the same value (flnite) as the mesh is reflned. 



We have run the same experiment for A(0) = 2.10947 and A(0) = 10. For both values the 
condition (5.12) is verifled, and therefore in view of Lemma 5.11 we expect the solution / to 



become discontinuous. To assess the ability of the DG methods to capture the singularity (or 
change in smoothness of /) we have measured the time evolution of the discrete and L^- 
norms of Vfh and study how they are affected under mesh reflnement and by an increment 
in the polynomial degree {k = 1,2,3,4). The results for A(0) = 2.10947 are depicted in 
Figure [23| 
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Figure 24. Plane Diode: Solution of the VP system for non-smooth solution 
test case at different times using mesh 60 x 60 and k = 2. 



Observe that the behaviour under mesh refinement of the measured quantities is very 



different for A(0) = (Figure [22]) and A(0) = 2.10947 (Figure [23]). From the graphics, for 
the former case we deduce convergence, indicating that the solution is indeed continuous. For 
A(0) = 2.10947, the ||V//j||l2 increases as the mesh is refined, which seems to indicate that 
1 1 V//^ 1 1^2 diverges (and blows up in finite time). The same effect is observed by increasing 
the polynomial degree. 



In Figure 24 we represent the approximate solution in phase space (x^v) for different times 



for A(0) = 2.10947. 



The corresponding results and graphics for A(0) = 10 are given in Figure 25 
in this case, by refining the mesh the quantity ||V//j||l2 and || V//j||l' 



Notice that 



also increases although 

after some time it decreases again. This effect might be due to the fact that the singularity 
in the problem is very weak, and as the time evolves the full discretized DG scheme might 
be adding too much artificial diffusion to capture the singularity at all times (note that the 
explicit time discretization adds also numerical diffusion) . The issue of coupling the scheme 
with some conservative integrator is currently under investigation. 

For A(0) = 10 we plot in Figure 26 the approximate density in phase space, computed with 
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Figure 25. Plane Diode: Evolution of ||V//^||l2 (top), ||//i||l2 (middle) and 
||V//^||L°° (bottom) for different mesh size and different k. 



a mesh 60 x 60 and polynomial degree k = 2. To further assess the possible ability of the 
DG schemes to capture the discontinuity in fh we have represented in Figure [27| the profile 
of the solution. 



6. Conclusions 



We have studied the verification and validation of the high order DG methods introduced 
in [1] for approximating the one-dimensional Vlasov-Poisson system with periodic bound- 
ary conditions. We have proposed two possible modifications of the definition of the DG 
schemes, that allow for practical and efficient implementation. We have shown theoretically 
and demonstrated numerically, that with such modifications the resulting DG methods still 
preserve the total discrete mass and energy (this if /c > 2). We have also verified numerically 
that in the case of smooth solutions, the approximate distribution function and the electro- 
static field converge optimally in L^. We have discussed the time integration for the schemes, 
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(c) t = 0.2 (d) t = 0.3 

Figure 26. Plane Diode: Solution of the VP system for non-smooth solution 
test case at different times using mesh 60 x 60 and k = 2. 



demonstrating that there is no essential benefit in using a TVD RK integrator (rather than 
a standard high order RK) for the practical simulations in plasma physics. We have shown 
how the fully discretized methods can be efficiently implemented in parallel. Moreover the 
performance of the introduced DG methods is validated, with several benchmark problems 
of plasma physics such as linear and nonlinear Landau damping and two different tests on 
two stream instability. We have also discussed how the schemes could be adapted for approx- 
imating a Vlasov-Poisson boundary value problem. It is demonstrated that the schemes have 
some potential ability for capturing the discontinuity of the solution in this case, although 
some further tuning on the time integration seems to be required to reproduce the correct 
weak singularity. This is currently under research. 
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Appendix A. Basis functions 



A.l. Lagrange polynomials. In the following parts we introduces basis functions that we 
used in the actual implementation of the scheme. Let / := [0, 1] then in order to span the 
P^(/) we introduce the Lagrange basis function defined by 

(A.l) 4(r):= n ^-^^ rel,yn = l,...,k + l, 

l<m</c+l,m7^n 

where {ri, ...^rk+i} is the set of distinct nodal coordinates in /. Recall that the Lagrange 
polynomials satisfy 

(A.2) L{rm) = Snm Vn, 771 = 1 , . . . , + L 
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Now let G P^(/) X P^(/), then we can write it as a tensor product 

k+l 

(A. 3) fh{x,V,t)= ^ C^n,m{t)ln{x)lm{v), 

n^m=l 

where an,m{t) are interpolation coefficients of fh{x^v^t) at (x^,^rn) and 
(A.4) {(^1, ^i), (^1, V2), (xfc+i, Vk+i)}, 

is the set of distinct nodal coordinates in P. For more information about the other pos- 
sible basis functions and methods to evaluate mass and gradient matrices using Lagrange 
polynomials we refer to [40] . 



A. 2. Bernstein polynomials. Here we briefly introduce a basis function for P^(/) that 
has a useful properties and will be used to facilitate evaluation of the term containing flux 
of electric fleld, i.e. E}^f}^^ in (2.2). Consider we are interested to determine whether or not 
Eh G P^(/) changes sign in / := [0, 1]. A simple approach to show that is positive (or 
negative) over whole / is to express E^ using polynomials that are positive in /. To do so 
we use Bernstein polynomials, B^{x)^ which are positive over / and satisfy 

(A.5) = (1 - x)B'^-\x) + X Bt\{x) Vn = 0, K 

for Bl{x) = 1 and B^{x) = 0, Vn < or n > fc. 



We use two properties of the Bernstein in the actual implementation. First one is 
(A.6) 5^(x)>0 Vn = 0, xg[0, 1]. 

More precisely, the Bernstein polynomials are all positive over /. Let us express Eh{x) by 

k 

then Eh{x) is non-negative over / if > 0,Vm (and similarly non-positive if < 0). 
Another properties says that Eh{x) is located inside the convex hull produced by the set 
{/^m}, which implies 

min^^ < min Eh{x), 
(A.7) ^ 

mSixEh{x) < max;5^. 

This last properties helps us to obtain an estimate of min and max of E^. For more details 
on other properties of Bernstein method that we do not use here, we refer to [42]. 

As we discussed before we use Lagrange polynomials to span P^ and therefore for E^. Now 
the question is how one can convert coeflicients of expansion using Lagrange basis functions 
to Bernstein. The change of basis from Lagrange to Bernstein involves a matrix-vector 
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multiplication, where the matrix is the inverse of a Vandermonde matrix. We expand Eh{x) 
in Lagrange polynomial basis 



(A.8) 



k+l 

m=l 



and define the vector of coefficients in Lagrange and Bernstein by 
(A. 9) f3 := [/3o, , a := [^i, ak+if , 

then we have 

(A.IO) cx = V 'f3, 

where the Vandermonde matrix is defined by 

(A.ll) Km:=5^_i(rn) Vn,mG l,...,fc + L 

If we want to check positivity (similarly negativity) of expressed using Lagrange polyno- 
mials (given a), we have to compute, /3 = • cx and then we just check positivity of all 
entries of /3 (and similarly for negativity). 

Note that since 1/ is a (/c + 1) • (/c + 1) Vandermonde matrix, its inversion is cheap and done 
once for all, before the time marching process starts and it is saved for later use. Therefore 



In Figure [28| we depict Lagrange, Legendre and Bernstein polynomials for k = 2. Note the 
positivity of Bernstein polynomials on / = [0, 1]. 



Lagrange, k = 2 



Legendre, k <2 



Bernstein, k = 2 
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Figure 28. Lagrange, Legendre and Bernstein basis functions of degree 2. 
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